library(MASS)
library(xtable)
library(tseries)
library(lmtest)
library(TSA)
library(vars)
library(fracdiff)
library(strucchange)
library(tsDyn)
library(fUnitRoots)


Data<-read.csv("rareearth_vekasi_data.csv")
dim(Data)
attach(Data)
names(Data)


#Quota
Japan_and_US_YTD<-(Japan_ChinaYTDQ+US_China_QYTD)*.001
quota.difference<-Quota-Japan_and_US_YTD

quota<-ts(quota.difference/1000,frequency=12,start=c(2005))
seasonal.quota<-stl(quota,s.window="per")
plot(seasonal.quota)
#extract trend
seasonal<-seasonal.quota$time.series[,"trend"]

quota.difference.1<-Quota-Japan_ChinaYTDQ*.001

quota.1<-ts(quota.difference.1/1000,frequency=12,start=c(2005))
seasonal.quota.1<-stl(quota,s.window="per")
plot(seasonal.quota.1)
#extract trend
seasonal.1<-seasonal.quota$time.series[,"trend"]



#Price variables-Japan

plot(ts(Japan_World_V,start=2005,frequency=12), col="red",ylab=)
lines(ts(Japan_China_V,start=2005,frequency=12),col="blue")


World.Value.Japan<-Japan_World_Q/Japan_World_V
China.Value.Japan<-Japan_China_Q/Japan_China_V

plot(ts(China.Value.Japan,start=2005,frequency=12), col="red",ylab=)
lines(ts(World.Value.Japan,start=2005,frequency=12),col="blue")

japan.price.difference<-World.Value.Japan/China.Value.Japan
japan.price<-ts(japan.price.difference,start=2005,frequency=12)

plot(ts(japan.price.difference,start=2005,frequency=12))

#Price variables-US
World.Value.US<-US_World_Q/US_World_V
China.Value.US<-US_China_Q/US_China_V

plot(ts(US_China_V,start=2005,frequency=12), col="red",ylab=)
lines(ts(US_World_V,start=2005,frequency=12),col="blue")



plot(ts(China.Value.US,start=2005,frequency=12), col="red",ylab=)
lines(ts(World.Value.US,start=2005,frequency=12),col="blue")

us.price.difference<-World.Value.US/China.Value.US
us.price<-ts(us.price.difference,start=2005,frequency=12)

plot(ts(us.price.difference,start=2005,frequency=12))


######MODEL

#Variables

diverse.japan<-ts(Japan_HHI_Normal, start=2005, frequency=12)
diverse.us<-ts(US_HHI_Normal, start=2005, frequency=12)
japan.price<-ts(japan.price.difference,start=2005,frequency=12)
us.price<-ts(us.price.difference,start=2005,frequency=12)
conflict<-ts(Exportban,start=2005,frequency=12)
conflict.quick<-ts(Shortban,frequency=12,start=c(2005))
quota.plain<-ts(Quota,frequency=12,start=2005)
seasonal<-seasonal.quota$time.series[,"trend"]

#First, look at Japan
#check for stationarity with augmented Dickey-Fuller tests
adf.test(seasonal)
adf.test(diverse.japan)#perhaps stationary, p=.04228
adf.test(japan.price)
#not stationary (closer than before)

#first differences
seasonal.diff.japan<-diff(seasonal)
diverse.diff.japan<-diff(diverse.japan)
price.diff.japan<-diff(japan.price)
#check for stationarity with ADFs of first differences
adf.test(seasonal.diff.japan)#not stationary, p=.217
adf.test(diverse.diff.japan)#stationary
adf.test(price.diff.japan)#stationary

#Phillips-Perron tests
pp.test(diverse.japan)#stationary
pp.test(japan.price)#stationary
pp.test(seasonal)#not stationary
#and now on the first differences
pp.test(diverse.diff.japan)#stationary
pp.test(price.diff.japan)#stationary
pp.test(seasonal.diff.japan)#not stationary

urzaTest(diverse.japan,model="trend")#shows that there is a potential structural break in the trend at position 31 (July, 2007)
urzaTest(japan.price,model="trend")#shows that there is a potential structural break in the trend at position 85 (January, 2012)
urzaTest(seasonal,model="trend")#shows that there is a potential structural break in the trend at position 112 (April, 2014)

urkpssTest(diverse.japan, lags="short")
urkpssTest(diverse.japan, lags="long")
urkpssTest(diverse.japan, lags="nil")

urkpssTest(japan.price, lags="short")
urkpssTest(japan.price, lags="long")
urkpssTest(japan.price, lags="nil")


urkpssTest(seasonal, lags="short")
urkpssTest(seasonal, lags="long")
urkpssTest(seasonal, lags="nil")


#are they cointegrated? Test using Philips-Ouliaris conintegration test
po.test(cbind(diverse.japan,japan.price))#yes (barely, p=.05748)
po.test(cbind(diverse.japan,seasonal))#yes
po.test(cbind(seasonal, japan.price))#no (p=.15, PO=-1.2779)
#Cointegration justifies use of regression model.
#should use a vecm, not a var.

endo<-data.frame(diverse.japan,japan.price,seasonal)
exo1<-data.frame(conflict)
exo2<-data.frame(conflict.quick)


allendo1<-data.frame(diverse.japan,japan.price,seasonal,conflict)
allendo2<-data.frame(diverse.japan,japan.price,seasonal,conflict.quick)


#test data for cointegration
philips.ouliaris.test<-ca.po(endo, demean="trend")
summary(philips.ouliaris.test)#test statistic 105.1473 (far above critical values)

#Johansen test, make cointegrated object
coint1<-ca.jo(endo,type="eigen",ecdet="trend",K=2,spec="longrun")
summary(coint1)

#restricted model
model.ecm1<-cajorls(coint1, r=1,reg.number=1)
summary(model.ecm1$rlm)
model.ecm1
par(mfrow=c(2,2))
plot(model.ecm1$rlm)


#unrestricted model
model.ecm1.2<-cajools(coint1, reg.number=1)
summary(model.ecm1.2)


coint2<-ca.jo(endo,type="eigen",ecdet="trend",K=2,spec="transitory")
summary(coint2)#very similar to longrun specification

#add dummy variable for conflict
conflict.dummy<-as.matrix(conflict)
colnames(conflict.dummy)<-"confdummy"

coint3<-ca.jo(cbind(diverse.japan,japan.price,seasonal),dumvar=conflict.dummy,ecdet="trend",K=2,spec="longrun")
summary(coint3)
#r still 1

model.ecm2<-cajorls(coint3,r=1,reg.number=1)#reg.number=1 estimates only for the diverse variable
summary(model.ecm2$rlm)
par(mfrow=c(2,2))
plot(model.ecm2$rlm)


#Try with lag of 3 as robustness check
coint3.1<-ca.jo(cbind(diverse.japan,japan.price,seasonal),dumvar=conflict.dummy,ecdet="trend",K=3,spec="longrun")
summary(coint3.1)
#r still 1 with 90% confidence

model.ecm2.1<-cajorls(coint3.1,r=1,reg.number=1)#reg.number=1 estimates only for the diverse variable
summary(model.ecm2$rlm)
par(mfrow=c(2,2))
plot(model.ecm2$rlm)




#short conflict dummy
conflict.dummy.short<-as.matrix(Shortban)
colnames(conflict.dummy.short)<-"confdummy.short"
coint4<-ca.jo(cbind(diverse.japan,japan.price,seasonal),dumvar=conflict.dummy.short,ecdet="trend",K=2,spec="longrun")
summary(coint4)
#r still 1

model.ecm3<-cajorls(coint4,r=1,reg.number=1)#reg.number=1 estimates only for the diverse variable
summary(model.ecm3$rlm)#pulse dummy not signficant

par(mfrow=c(2,2))
plot(model.ecm3$rlm)



#Next, look at United States
#check for stationarity with augmented Dickey-Fuller tests
adf.test(seasonal)
adf.test(diverse.us)
adf.test(us.price)
#not stationary 

#first differences
seasonal.diff.us<-diff(seasonal)
diverse.diff.us<-diff(diverse.us)
price.diff.us<-diff(us.price)
#check for stationarity with ADFs of first differences
adf.test(seasonal.diff.us)#not stationary, p=.217
adf.test(diverse.diff.us)#stationary
adf.test(price.diff.us)#stationary

#Phillips-Perron tests
pp.test(diverse.us)#stationary
pp.test(us.price)#stationary
pp.test(seasonal)#not stationary
#and now on the first differences
pp.test(diverse.diff.us)#stationary
pp.test(price.diff.us)#stationary
pp.test(seasonal.diff.us)#not stationary

urzaTest(diverse.us,model="trend")#shows that there is a potential structural break in the trend at position 87
urzaTest(us.price,model="trend")#shows that there is a potential structural break in the trend at position 73
urzaTest(seasonal,model="trend")#shows that there is a potential structural break in the trend at position 76

urkpssTest(diverse.us, lags="short")
urkpssTest(diverse.us, lags="long")
urkpssTest(diverse.us, lags="nil")

urkpssTest(us.price, lags="short")
urkpssTest(us.price, lags="long")
urkpssTest(us.price, lags="nil")


urkpssTest(seasonal, lags="short")
urkpssTest(seasonal, lags="long")
urkpssTest(seasonal, lags="nil")


#are they cointegrated? Test using Philips-Ouliaris conintegration test
po.test(cbind(diverse.us,us.price))#yes
po.test(cbind(diverse.us,seasonal))#yes
po.test(cbind(seasonal, us.price))#borderline (p=.08165)
#Cointegration justifies use of regression model.
#should use a vecm, not a var.

endo.us<-data.frame(diverse.us,us.price,seasonal)
exo1<-data.frame(conflict)
exo2<-data.frame(conflict.quick)


allendo1.us<-data.frame(diverse.us,us.price,seasonal,conflict)
allendo2.us<-data.frame(diverse.us,us.price,seasonal,conflict.quick)


#test data for cointegration
philips.ouliaris.test.us<-ca.po(endo.us, demean="trend")
summary(philips.ouliaris.test.us)#test statistic 73.1316 (above critical values)

#Johansen test, make cointegrated object
coint1.us<-ca.jo(endo.us,type="eigen",ecdet="trend",K=2,spec="longrun")
summary(coint1.us)


#restricted model
model.ecm1.us<-cajorls(coint1.us, r=1,reg.number=1)
summary(model.ecm1.us$rlm)
model.ecm1.us
par(mfrow=c(2,2))
plot(model.ecm1.us$rlm)


#unrestricted model
model.ecm1.2.us<-cajools(coint1.us, reg.number=1)
summary(model.ecm1.2.us)


coint2.us<-ca.jo(endo.us,type="eigen",ecdet="trend",K=2,spec="transitory")
summary(coint2.us)#very similar to longrun specification

#add dummy variable for conflict
conflict.dummy<-as.matrix(conflict)
colnames(conflict.dummy)<-"confdummy"

coint3.us<-ca.jo(cbind(diverse.us,us.price,seasonal),dumvar=conflict.dummy,ecdet="trend",K=2,spec="longrun")
summary(coint3.us)
#r still 1, now at 99% confidence

model.ecm2.us<-cajorls(coint3.us,r=1,reg.number=1)#reg.number=1 estimates only for the diverse variable
summary(model.ecm2.us$rlm)
par(mfrow=c(2,2))
plot(model.ecm2.us$rlm)


#short conflict dummy
conflict.dummy.short<-as.matrix(Shortban)
colnames(conflict.dummy.short)<-"confdummy.short"
coint4.us<-ca.jo(cbind(diverse.us,us.price,seasonal),dumvar=conflict.dummy.short,ecdet="trend",K=2,spec="longrun")
summary(coint4.us)
#r still 1

model.ecm3.us<-cajorls(coint4.us,r=1,reg.number=1)#reg.number=1 estimates only for the diverse variable
summary(model.ecm3.us$rlm)#pulse dummy not signficant

par(mfrow=c(2,2))
plot(model.ecm3.us$rlm)




